!Dynamique moléculaire d'annihilation 1D avec condition initiale bimodale: \ln(N/N_now) \sim - \xi \ln(t) 
!But: transition entre régime en 1/\sqrt(t) [temps court - corrélations] et 1/t [temps longs - choas moléculaire]
!Nouvel essai: 20-09-2003

module fonctions
implicit none
integer, parameter :: N = 10000 !nombre initial de particules
integer, parameter :: avg = 3 !nombre de runs pour moyenne
real, parameter :: intervalle = 1.e-4 !intervalle pour plot log-log
character*(*), parameter :: filename = 'MDbimodal1d.txt' !nom du fichier de sortie
double precision, parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164
contains

FUNCTION ran2(idum)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL :: ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,&
          &IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER :: idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if(idum.le.0)then
	idum=max(-idum,1)
	idum2=idum
	do j=NTAB+8,1,-1
		k=idum/IQ1
		idum=IA1*(idum-k*IQ1)-k*IR1
		if(idum.lt.0)idum=idum+IM1
		if(j.le.NTAB)iv(j)=idum
	enddo
	iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if(idum.lt.0)idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if(idum2.lt.0)idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
end function ran2

!----------------------------------------------------------------------------------

SUBROUTINE hpsort(n,m,ra) !size: ra(1:n,1:m), sort according to ra(:,1)
      INTEGER n,m
      DOUBLE PRECISION ra(n,m)
      INTEGER i,ir,j,l
      DOUBLE PRECISION rra(m)
      if (n.lt.2) return
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
          l=l-1
          rra=ra(l,:)
        else
          rra=ra(ir,:)
          ra(ir,:)=ra(1,:)
          ir=ir-1
          if(ir.eq.1)then
            ra(1,:)=rra
            return
          endif
        endif
        i=l
        j=l+l
20      if(j.le.ir)then
          if(j.lt.ir)then
            if(ra(j,1).lt.ra(j+1,1))j=j+1
          endif
          if(rra(1).lt.ra(j,1))then
            ra(i,:)=ra(j,:)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        goto 20
        endif
        ra(i,:)=rra
      goto 10
END SUBROUTINE hpsort

end module fonctions

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------

program MDbimodal
use fonctions
implicit none
double precision, allocatable, dimension(:,:) :: particules !particules(x,1) = position, particules(x,2) = vitesse
double precision, dimension(1:N,1:4) :: d !1: distance relative, 2: vitesses relatives, 3-4: voisins gauche-droite
double precision, dimension(1:4) :: tmp
double precision :: time,distmin,told
real, dimension(:,:), allocatable :: result
integer :: i,j,seed,Nnow,run,sizemin,pos

seed = -14731

i=nint(1.2/intervalle)
allocate(result(1:i,1:2))
sizemin = 100000000
do run=1,avg !boucle pour moyenne sur runs

allocate(particules(1:N,1:2))
particules = 0.d0
do i=1,N !initialisation sur la droite [0,1]
	particules(i,1) = dble(ran2(seed))
	if (ran2(seed) > 0.5) then
		particules(i,2) = 1.d0 ! ->
	else
		particules(i,2) = -1.d0 ! <-
	end if
end do

call hpsort(N,2,particules) !met dans l'ordre croissant les particules

do i=1,N-1 !espace des distances relatives. Vitesses: 0: -> <- 1: -> -> -1: <- ->
	d(i,1) = particules(i+1,1) - particules(i,1)
	if (nint(particules(i,2)) == 1) then
		if (nint(particules(i+1,2)) == -1) then
			d(i,2) = 0.d0
		else
			d(i,2) = 1.d0
		end if
	else
		if (nint(particules(i+1,2)) == -1) then
			d(i,2) = 1.d0
		else
			d(i,2) = -1.d0
		end if
	end if
end do
d(N,1) = (1.d0 - particules(N,1)) + particules(1,1)
if (nint(particules(N,2)) == 1) then
	if (nint(particules(1,2)) == -1) then
		d(N,2) = 0.d0
	else
		d(N,2) = 1.d0
	end if
else
	if (nint(particules(1,2)) == -1) then
		d(N,2) = 1.d0
	else
		d(N,2) = -1.d0
	end if
end if
do i=2,N-1 !voisins
	d(i,3) = dble(i)-1.d0
	d(i,4) = dble(i)+1.d0
end do
d(1,3) = dble(N)
d(1,4) = dble(2)
d(N,3) = dble(N-1)
d(N,4) = dble(1)
deallocate(particules)

pos = 0
Nnow=N
told = 1.e-7
time = 1.d0
open(unit=1,file=filename)
close(1,status='DELETE') 
write(*,*)
write(*,'("Run:",i3,"/",i3)') run,avg
onerun: do !décimation de 1 run

	if (all(nint(d(1:Nnow,2)) == 1)) exit onerun !toutes les particules vont dans la même direction

	distmin=10.d0 !trouve la plus petite distance avec collision
	findfirst: do i=1,Nnow
		if (d(i,1) < distmin) then
			if (nint(d(i,2))==0) then
				distmin = d(i,1)
				j=i
			end if
		end if
	end do findfirst

	if (j<Nnow .and. j>1) then !nouvelle distance, nouveau signe, nouveaux voisins de j
		tmp(1) = d(j-1,1) + d(j+1,1) !nouvelle distance
		if (nint(d(j-1,2)) == 1 .and. nint(d(j+1,2)) == 1) then !nouveau signe
			tmp(2) = 0.d0
		else
			tmp(2) = - d(j-1,2)*d(j+1,2)
		end if
		tmp(3) = d(j-1,3) !nouveaux voisins de j
		tmp(4) = d(j+1,4)
		where (nint(d(:,4)) == j-1) d(:,4) = dble(j) !nouveaux voisins des voisins de j
		where (nint(d(:,3)) == j+1) d(:,3) = dble(j)
	else
		if (j==1) then
			tmp(1) = d(Nnow,1) + d(2,1) !nouvelle distance
			if (nint(d(Nnow,2)) == 1 .and. nint(d(2,2)) == 1) then !nouveau signe
				tmp(2) = 0.d0
			else
				tmp(2) = - d(Nnow,2)*d(2,2)
			end if
			tmp(3) = d(Nnow,3) !nouveaux voisins de j
			tmp(4) = d(1,4)
			where (nint(d(:,4)) == Nnow) d(:,4) = 1.d0 !nouveaux voisins des voisins de j
			where (nint(d(:,3)) == 2) d(:,3) = 1.d0
		else
			tmp(1) = d(Nnow-1,1) + d(1,1)
			if (nint(d(Nnow-1,2)) == 1 .and. nint(d(1,2)) == 1) then !nouveau signe
				tmp(2) = 0.d0
			else
				tmp(2) = - d(Nnow-1,2)*d(1,2)
			end if
			tmp(3) = d(Nnow-1,3) !nouveaux voisins de j
			tmp(4) = d(1,4)
			where (nint(d(:,4)) == Nnow-1) d(:,4) = dble(Nnow) !nouveaux voisins des voisins de j
			where (nint(d(:,3)) == 1) d(:,3) = dble(Nnow)
		end if
	end if

	time = time + d(j,1) !corrige toutes les autres distances
	where (nint(d(1:Nnow,2)) == 0) d(1:Nnow,1) = d(1:Nnow,1) - d(j,1)
	where (nint(d(1:Nnow,2)) == -1) d(1:Nnow,1) = d(1:Nnow,1) + d(j,1)

	d(j,:) = tmp !nouvel intervalle et autres

	if (j<Nnow .and. j>1) then !Met les deux distances j-1 et j+1 la fin du tableau
		d(j-1,:) = d(j,:)
		d(j:Nnow-2,:) = d(j+2:Nnow,:)
	else
		if (j==1) then
			d(2:Nnow-2,:) = d(3:Nnow-1,:)
		else
			d(1:Nnow-3,:) = d(2:Nnow-2,:)
			d(Nnow-2,:) = d(Nnow,:)
		end if
	end if
	d(Nnow-1:Nnow,:) = 100000.d0

	Nnow = Nnow-2 !actualise le nombre de particules
	if (Nnow <= 6) exit onerun

	if (modulo(real((N-Nnow)*100)/real(N),1.)==0.) write(*,'(i3,"%")', ADVANCE='NO') nint(real((N-Nnow)*100)/real(N))
	if (dlog10(1.d0+time-told) >= intervalle) then
		told = time
		pos = pos+1
		result(pos,1) = result(pos,1) + real(time)
		result(pos,2) = result(pos,2) + real(Nnow)/real(N)
	end if

end do onerun
if (pos < sizemin) sizemin = pos
end do

result = result/real(avg)
open(unit=1,file=filename)
do i=1,sizemin
	write(1,*) result(i,1), result(i,2)
end do
close(1)
deallocate(result)
write(*,*) "End."
read(*,*)
end program MDbimodal